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ABSTRACT 

An algorithm is developed and tested for the problem posed by photometric obser- 
vations of the bulge of the Milky Way. The latter subtends a non-trivial solid angle 
on the sky, and we show that this permits inversion of the projected brightness dis- 
tribution under the assumption that the bulge has three orthogonal mirror planes of 
specified orientation. A serious error in the assumed orientation of the mirror planes 
should be detectable. 
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1 INTRODUCTION 

In the last few years evidence has been accumulating that 
the galactic bulge is barred (e.g. Blitz 1993; Gerhard 1995). 
The distribution of near-infrared light from the direction 
of the galactic centre provides the most direct evidence for 
non-axisymmetry of the bulge since at A ~ 2 /tin starlight 
is the principal contributor to the overall surface brightness 
at (|l| < 20°, \b\ < 8°), and at such wavelengths our view of 
this region is only moderately obscured by dust. 

Blitz & Spergel (1991) first detected tell-tale signs of 
the bulge's triaxiality in the 2 /tin data of Matsumoto et 
al (1982), but these old data were far too noisy to permit 
detailed model-building. The near-infrared brightness dis- 
tribution determined by the DIRBE experiment aboard the 
COBE satellite (Weiland et al 1994) are of very much higher 
quality and several groups are currently basing models of 
the bulge on these data. This activity raises the following 
question: how uniquely does the projected brightness distri- 
bution of the bulge determine the three-dimensional lumi- 
nosity density in the bulge? We address this question. We 
are primarily interested in the feasibility, in principle, of the 
inversion of ideal data, rather than the practical question of 
whether any given body of data is of sufficient quality to be 
usefully inverted. 

Since the bulge's projected brightness distribution is a 
function of the two coordinates of the sky, while the lumi- 
nosity density of a bar is a function of all three spatial coor- 
dinates, it would be natural to assume that the latter cannot 
be unambiguously determined from the former. We show, 
however, that it is possible in many cases to recover ac- 
curately the three-dimensional luminosity distribution from 
the projected density provided the bulge is known to have 
three mutually orthogonal mirror planes of known orienta- 
tion. Thus the ambiguity inherent in determination of the 
luminosity density of the bulge can be largely reduced to 
that involved in choosing the orientation of the bulge's mir- 
ror planes. Moreover, it turns out that the algorithm we 
develop can offer valuable guidance as to the orientation of 
these planes. 



In Section 2 we formulate a Richardson-Lucy (Richard- 
son 1972; Lucy 1974) deprojection algorithm for the recovery 
of bulge models from surface-brightness data, and explain 
how the assumption that the bulge is eight-fold symmetric 
about given axes enables a three-dimensional model to be 
inferred from two-dimensional data. In Section 3 we test 
the performance of this algorithm on pseudo-data. Section 
4 sums up. 

2 A RICHARDSON LUCY ALGORITHM 
2.1 Motivation 

As Contopoulos (1954) already showed, the elliptical image 
of a distant galaxy can be interpreted as the projection of an 
infinite number of different elliptical bodies (see also Stark 
1977). In particular, the image of any body whose isodensity 
surfaces are self-similar ellipsoids can be interpreted as the 
image of an axisymmetric body. It is helpful to understand 
why the non-negligible angular size of the Milky Way's bulge 
may allow less freedom in the interpretation of its photom- 
etry, and in particular betray the bulge's triaxiality. 

From photometry alone, one cannot tell from what dis- 
tance one receives any given contribution to the brightness 
observed at some point on the sky. Consequently, one can 
deform a galaxy into a different object with the same image 
by moving stars along their lines of sight. (In general the lu- 
minosities of the stars will have to be adjusted as they move 
nearer to or further from the observer.) When a galaxy is 
viewed from afar, all lines of sight are parallel and such re- 
organizations of the galaxy can easily carry it between two 
states of symmetry, such as ellipsoidal brightness distribu- 
tions. 

When a galaxy subtends a non-trivial solid angle, it is 
not clear that the system can be deformed between states of 
symmetry without changing the projected image. In partic- 
ular, if we wish a deformation to be compatible with the sys- 
tem retaining mirror symmetry in three orthogonal planes, 
motion of any star along its line of sight will have to be 
accompanied by the motion of its mirror-image stars in di- 
rections that generally do not run parallel to their lines of 
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sight. Thus such deformations of the three-dimensional fig- 
ure of the galaxy would be accompanied by deformation of 
the image. The mathematical scheme below exploits this 
connection to infer a mirror-symmetric distribution from a 
projected image. 

2.2 Mathematical formulation 

We employ coordinate axes which coincide with the principal 
axes of the bulge. Then associated with any point x ^ o, 
there are eight distinct points x t , i = 1, . . . , 8, symmetrically 
placed with respect to the principal axes, such that v(xi) = 
v(x), where v[x) is the luminosity density due to stars at 
x. We refer to the the "symmetric companions" of x. 

The assumed octant-symmetry of the bulge may be ex- 
pressed mathematically by 



v{x) 



where 



d 3 x' 8g(x, x')v(x'), 



8 s (x,x') = 6(| zi | — xj)6(|x 2 | — x' 2 )8(\x3 \ — x' 3 ). 



(1) 



(2) 



Notice that 8g vanishes unless its second argument lies in 
the positive octant, x t > 0. 

The surface brightness in the direction from the Sun J7 

is 



/(J?) = / dsu(x(s)), 
lo 



(3) 



where x(s) = Xq + sfl is the point distance s from the Sun 
in the direction fl. Substituting from (1) into (3), we have 



J(J7) = J d 3 iV(i') J As8g(x(s),x') 
= / & 3 x'v(x')K(fl\x'), 

J -foctant 

where 

l> CO 

K(fi\x') = / ds8g(x(s),x'). 



(4a) 



(4b) 



The kernel K does not satisfy the normalization condi- 
tion J d 2 fi K = 1 required of a probability distribution. In 
fact 



& 2 flK{n\x') = / d 2 J7 / ds8g(x(s),x') 



E 



s 2 (x' t ) - l 2 (x') 



for x' not 
in +octant, 

otherwise. 



Here x' t is the position of the i th symmetric companion of 
x' and s(x') = \x' — xq\ is the distance from the Sun to x' . 
In view of (4) we define 

K(fi\x') = l 2 {x')K{fl\x') and v{x') = ^fy ■ (5) 

s (x ) 

K(fl\x) obviously vanishes unless the direction J7 points to 
one of the x t . When fl does so point, K is ^-function like. 



Integrating K over small a solid angle AJ7 8 around fl t , we 
have 



& 2 riK(n\x') = s (x ) I A 2 nK(n\x' 



s 2 (x') J d 2 J7 J &s8g(x(s),x'^ 
_ 2 , f f d 2 J7s 2 ds , 



s (x ) 
s*(x>) 



s 2 {x l ) 



-8g (x(s), x') 



£Li/* 2 K-) * 2 (*l) 



It follows that 



_ 2 . ^ 8(n- n 
K{n\x) = I s ( x > 



3 2( Xi ) for x in + octant (6) 

i = l 



v = otherwise. 
In conclusion, our Richardson-Lucy scheme for the in- 
version of equation (3) is 



(7) 



I r (fl)= I d xv r (x)K(fl\x) 
d 3 x v r (x)K( fi\x) 
ds v T (x(s)^ , 



as expected, and 

v r+1 (x) = v r (x) / d 2 riK(n\x) 



I r (fl) 



v r (x) j d 2 J7 s 2 (x) y ^ 



g (8) 

8{n - n t ) i{fi) v ; 



s 2 (xi) i r (ny 



Or, simply 



V r +l(x) = V r (x)s 2 (x) T 



I{fl t ) 1 



{ni)s 2 {xiy 



(9) 



The convergence properties of any R-L scheme may be 
improved by modifying the scheme so that the 'data' are not 
I(fl) but c/(J7)/(J7), where g(fi) is an arbitrary positive 
function (Lucy 1994; Dehnen 1995). The iterates I r pro- 
duced by the modified scheme are such that gl In gl r rather 
than Ilnl r is maximized as r — > oo. Multiplying both sides 
of equation (4a) by g, it is clear that the modified scheme 
has kernel gK , and it is easy on to show that for this kernel 
(9) becomes 



(10) 



A ^ IjOj) g{n t ) gjOj) 

Vr+1 (x) = Vr {x) L Tm^W) -L 

i=i i=i 

In the experiments described below, the minimization of the 
residuals on the sky was usefully accelerated by taking g to 
be of the form g = I -13 with /? > 0; five iterations with /? = 
0.2 were followed by four iterations with fi = 0.8. The early 
iterations emphasize the residuals near the centre, while the 
later ones concentrate on diminishing the residuals further 
out. 
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3 TESTS OF THE ALGORITHM 

We have tested the effectiveness of the deprojection algo- 
rithm given by (10) as foiiows. For each point of an 80 x 80 
grid of points in the (I, b) plane we calculate the line-of-sight 
integral I = J ds v through a simple ellipsoidal density dis- 
tribution v[x). For most of our experiments v was of the 
form 



v(x) = (XeVeVe) exp 



2 i 2 

r c + x 



1/2 



(11) 



where r c = 175 pc is comparable with the spatial resolution 
of our calculation (see below). The line of sight from the 
Sun to x = is inclined at angle do to the (x, y) plane, and 
the x-axis is inclined at angle <j>o to the projection of the 
Sun-centre line onto the (x, j/)-principal plane. We adopt 
\xq\ = 8 kpc and integrate from s = 2 — 14kpc down each 
line of sight. Subsequently, the value of I at an arbitrary 
point (I, b) is inferred by bilinear interpolation between the 
values taken by In I on this sky-grid. 

Our model density distribution is defined by the value of 
v T at each point on a Cartesian grid of 25 x 25 x 20 points that 
covers the positive octant uniformly out to x max = 4.2 kpc, 
?/max = 3.2 kpc, z ml! = 2 kpc. We refer to the region \x\ < 
£max, \y\ < 2/max, \z\ < 2 m ax as the model's "bounding box." 
The model may have a different orientation from that of the 
true distribution v\ the Sun-centre line is inclined at angle 
9 to the model's equatorial plane, and the projection of the 
Sun-centre line into that plane is inclined at angle <f> with 
respect to the model's x-axis. 

For each point (I, b) of the sky-grid we project v by 
integrating along the section of the line of sight that lies 
within the bounding box. When performing J dsiv, values 
of v T are obtained by tri-linear interpolation of In v T on the 
three-dimensional grid. 

The sky-grid is chosen to be sufficiently large that it 
covers all points at which the projected density is non-zero. 
Since the bounding box projects to an irregular shape on the 
sky - its silhouette; see Fig. 1 - the sky-grid inevitably con- 
tains many lines of sight which do not intersect the bounding 
box. Along such lines of sight the projected density is set to 
an arbitrary small value. 

At points on the sky which lie near the boundary of the 
box's silhouette, the projected density is artificially small 
because only a small portion of the line of sight intersects 
the bounding box. Allowance must be made for this when 
comparing with either real data or our pseudo-data, which 
are obtained by integrating over the most significant 12 kpc 
of every line of sight. We make this allowance by projecting 
our initial model, which is given by an analytic function 
such as (11), twice, once along 12 kpc of each line of sight 
and once through the bounding box only. The ratio l\ and 
Iq between these two estimates is recorded and subsequently 
used to estimate I r through 



I - —I h 

1 



(12) 



where I r is the projection of v T through the bounding box. 
The best way of seeing that our R-L algorithm will reduce 
the residuals in I/I r with I r obtained in this way is to ob- 




Fig. 1 A trial at general orientation: 9q = = 7°, 0o = = 18°. 
After nine iterations an initially spherical model (dotted curves) 
has been deformed into a reasonable approximation (dashed 
curves) to the true surface-density contours (full curves). The 
contours of Iq crowd at the edge of the silhouette of the bound- 
ing box because outside this silhouette Iq is identical with Iq. 



serve that one can write 
I (lo/Io)l 



(13) 



T J h 

That is, the factor in brackets on the top of the right-hand 
side of (13) is used to estimate the fraction of the data I that 
arises in the bounding box, and then the R-L algorithm can 
be considered to be applied to projection from within the 
bounding box only. 

Fig. 1 shows the result of an experiment in which 
the orientations of the trial and true distributions were 
general but identical: do = = 7°, <j>o = <j> = 18°. 
[These values are motivated by the work of Burton & Liszt 
(1978), who find > 7°, and of Binney et al. (1991), who 
find <j> ~ 18°.] The true distribution had principal axes 
(i e , ?/e, z e ) = (0.5, 0.3, 0.15) kpc, while the first trial distribu- 
tion was spherical. The full and dotted curves show, respec- 
tively, the true and first-guess isophotes. The dashed curves 
show the isophotes of Ig. Near the centre these isophotes 
closely follow the true isophotes, but further out material 
deviations from the true isophotes are apparent, especially 
above the origin. The silhouette of the bounding box is de- 
lineated by the sharp crowding of the isophotes of Ig where 
the transition is made from the region in which they have 
been adjusted by the R-L algorithm to that in which they 
coincide with the original (dotted) isophotes. 

This sharp interface between the regions of modified and 
unmodified surface brightness arises because the initial guess 
(of a spherical distribution) was poor, with the result that 
the ratio P / Iq used in equation (12) is sub-optimal. In such 
a case it would be immediately apparent from the data that 
the true distribution is somewhat flattened, and one would 
adopt a flattened distribution for vq. Fig. 2 is the same as 
Fig. 1 but with vo given by equation (11) with x e = y e = 
0.4 kpc, z e = 0.15 kpc. The silhouette of the bounding box is 
now much less clear and the significant deviations between 
dashed and full isophotes are now confined to the upper left 
corner of the image, where the effects of truncation by the 
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Fig. 2 The same as Fig. 1 except that the initial distribution 
was flattened rather than spherical [vq was given by (11) with 
x e = y e = 0.4kpc, 2: e = 0.15kpc]. In consequence the correc- 
tion factors Iq/Iq are more accurate and the silhouette of the 
bounding box is less apparent. 
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Fig. 3 The density distributions that underlie the trial of Fig. 2. 
The full curves show the intersections of the true isodensity sur- 
faces with the (x, y) plane, while the dotted and the dashed curves 
show the same curves for vq and ^9 . 

bounding box are most severe. 

Fig. 3 shows, for the inversion of Fig. 2, the contours in 
which isodensity surfaces of the true model (full contours), 
the initial-guess model vq (dotted contours) and the iterated 
model i/g cut the (x,y) plane. The innermost half of the 
dashed contours follow the full contours quite accurately. 
Thus the algorithm has correctly deduced that the bulge is 
bar-shaped. Further out, the dashed contours follow the full 
contours for x < 3 kpc but then turn upwards towards the 
boundary. This effect arises from the need indicated above, 
to increase the model's surface brightness in the upper left 
of Fig. 2. In particular, if the iterations are started from 
an over-elongated distribution, such as (11) with axes x e = 
0.6 kpc, y e = 0.3 kpc, z e = 0.15 kpc, truncation by the box 
has the opposite effect, and the dashed curves in a plot such 
as Fig. 3 lie below the full curves at large x. Fig. 4 illustrates 
this point by showing the same physical region as is shown 
in Fig. 3 taken from a computation which employed a larger 
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Fig. 4 The same as Fig. 3 except that the bounding box has been 
extended to 6kpc, 4.6kpc and 2.9kpc in x, y and z, respectively, 
with a corresponding increase in the number of grid points within 
the box. The same physical region is shown. 

bounding box. The spurious effects of the bounding box 
have now moved outwards. 

Fig. 5 shows the result of a trial in which the prin- 
cipal axes of the density distribution that generated the 
data were rotated in the (x,y) plane with respect to the 
assumed symmetry axes of the deprojected model. Specif- 
ically, 8 = 9 = 7°, 4> = 30°, 4> = 18°. Fig. 6 shows that 
the algorithm now has difficulty diminishing the residuals 
beyond a certain point. Moreover, in Fig. 5 the residuals 
form a clear pattern; the predicted and true isophotes cross 
one another three times in each revolution of the projected 
centre. Clearly, one could search for the true orientation of 
the bar by finding the orientation at which the residuals on 
the sky are smallest and have no clear angular pattern. 

From Figs 1-5 we conclude that for #0 =0 = 7°, our 
R-L algorithm is able to infer the true brightness distribu- 
tion of a triaxial bulge to good accuracy. Spurious effects 
generated by the bounding box are readily diagnosed and 
can be used to guide one's initial guess towards one that 
lies sufficiently close to the truth to give good final results. 
Alternatively, they can be shifted out of the region of astro- 
physical interest by using a larger bounding box. When the 
assumed symmetry axes are inclined with respect to the true 
ones, the residuals on the sky develop a pattern which sug- 
gests that the model needs to be re-oriented. Similar results 
are obtained for other values of 9 and <f>, provided neither 
angle is chosen such that the Sun-centre line lies almost in 
a principal plane. 

Fig. 7, which shows the result of a trial with #0 = = 
and <j>o = <j> = 18° , shows that the scheme performs less well 
when the Sun-centre line lies close to the equatorial plane. 
Fig. 6 shows that the residuals in the (I, b) plane decline 
at exactly the same rate as in the test shown in Fig. 2, 
yet along a line in Fig. 7 that is inclined at ~ 18° to the 
x-axis, the dashed curves, which mark the intersection of 
the plane with the equidensity surfaces of 1/9 , turn abruptly 
downwards. As the inclination <f> of the model to the Sun- 
centre line is varied, the inclination of the line on which the 
dashed curves have shoulders follows <f>. The sharpness of 
the shoulders gradually reduces as 9 is increased. It seems 
that the algorithm correctly infers that v is elongated along 
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Fig. 5 A trial in which 9 = 9 = 7°, <f> = 30°, <f> = 18°. 
The residuals tend to an rms value of ln(/// r ) ~ 0.4 and 
a characteristic mismatch of projected and true isophotes is 
apparent. The true isophotes are shown by full curves and 
the model ones by dashed curves. The isophotes of the first 
(axisymmetric) model are shown dotted. This model had 

(x e ,ye,Ze) = (0.4, 0.4, 0.1 5) kpc. 
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Fig. 6 RMS residuals of ln(// I r ) for the four deprojections shown 
in other figures. Triangles correspond to Fig. 1, squares to Figs 2 
and 3, pentagons to Fig. 5 and hexagons to Fig. 7. 

the x-axis, but at small 9 the correct projected density can 
be produced with very non-elliptical isodensity surfaces, and 
for some reason the algorithm picks out such non-elliptical 
surfaces. 

When the model is viewed from within a principal plane, 
the images of the symmetric companions that lie above and 
below the plane no longer provide independent information 
on the model's light distribution - the assumed symmetry 
of the model and data already suffice to ensure that if model 
and data agree at b > 0, then they must also agree at b < 
0. Moreover, when 9 = 0, any motion of a source that 
leaves fixed on the sky the images of both the source and 
those of its symmetric companions that lie above the plane, 
will also leave fixed on the sky its symmetric companions 
below the plane. By contrast when 9 ^ 0, the symmetric 
companions below the plane will move on the sky even when 
those above the plane do not, and the amount by which they 
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Fig. 7 As for Fig. 3 except that the Sun-centre line lay in the 
equatorial plane (9q = 9 = 0, </>o = 4> = 18°). The true and 
model isodensity surfaces intersect the equatorial plane in the 
full and dashed curves, respectively. The isodensity surfaces of 
the initial model generated the dotted curves. The dashed curves 
have unphysical shoulders that lie along a line that is inclined to 
the ar-axis by an angle roughly equal to </>. 



will move will increase continuously with 9. Consequently, 
the information content of the data increases with 9. 

As we have remarked, we cannot expect the deprojec- 
tion to be unique even in the case 9 ^ 0, so it is not surpris- 
ing that in the degenerate case 9 = the algorithm picks a 
solution that differs materially from the source of the data. 
Fortunately, the spurious structure in the solution is readily 
identified because it is associated with the Sun-centre line: 
the difference between the original density and the recovered 
density consists largely of a cone of negative density around 
the x-axis and a ridge of positive density around the sun- 
centre line. This structure is reminiscent of axisymmetric 
konus densities (Gerhard & Binney 1995), which are inti- 
mately related to the non-uniqueness of the axisymmetric 
deprojection problem. 

One can ensure the data are fitted by more nearly el- 
liptical density distributions v r by restarting the algorithm 
from an elongated rather than a spherical distribution af- 
ter one has discovered that the data require elongated con- 
tours in the (x,y) plane. Alternatively, sharp features in 
the model can be discouraged by either filtering the model 
at each iteration or by biasing the Lucy iterations against 
them (Lucy 1994). 



3.1 Noise 

Since real data are always noisy, it is important to under- 
stand how robustly the algorithm copes with noise. Fig. 8 
shows the results of deprojecting two noisy data sets. In the 
left panel the problem was identical to that of Figs 2 and 
3 except that Gaussian noise was added to the data with 
dispersion a = 0.001/(0,0). This noise dominates the data 
for I > 20°, but rapidly becomes unimportant further in. 
The left panel in Fig. 8 shows that the density inside 3 kpc 
is accurately recovered, and outside this radius the model 
rapidly becomes worthless. The right panel in Fig. 8 shows a 
similar experiment in which the true luminosity density had 
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Fig. 8 Deprojections of noisy data. The left panel is as for Fig. 3 except that Gaussian noise has been added to the data at a level that 
dominates the surface brightness beyond / ~ 20° — the dispersion was 0.001/(0,0). The right panel shows the result of deprojecting a 
much shallower brightness distribution — the scale lengths in (11) were three times those used for the other figures. Again Gaussian noise 
[<r = 0.07/(0, 0)] dominates the signal outside / ~ 20°, but is now also important further in. These models were obtained after just four 
iterations with f3 = 0.2. 



much larger scale lengths. Specifically, in this case the pa- 
rameters (i e , y e , z e ) of equation (11) were (1.5, 0.9, 0.45) kpc. 
The noise level was again chosen such that noise and data 
were equal at I ~ 20°. In this case noise is important even at 
1 — 3 kpc. Nonetheless, the right panel in Fig. 8 shows that 
after four iterations the algorithm has correctly diagnosed 
the elongation of the bar. In a practical case, this very noisy 
inversion would inspire one to repeat the analysis, this time 
starting from an appropriately elongated model distribution. 
After one or at most two iterations one's model would then 
provide a good fit to the data while still not being strongly 
contaminated by noise. 

4 CONCLUSIONS 

We have presented a Richardson-Lucy algorithm for the re- 
covery of a mirror-symmetric density distribution from its 
image when viewed from a distance similar to our distance 
from the bulge of the Milky Way. 

Practical application of this algorithm is complicated 
by the necessity of somehow artificially confining the three- 
dimensional model being recovered. We place the model 
in a box and use the simple rule (12), which depends on 
the initially guessed distribution i>o, to take into account 
emission from outside the box. Since the quality of vq soon 
manifests itself in the predicted distribution on the sky, it is 
in practice straightforward to iterate towards a good choice 
for vo- 

When the assumed orientation of the principal axes dif- 
fers significantly from the true orientation, the algorithm 
correctly recovers the elongation of the isodensity surfaces, 
but the residuals on the sky remain relatively large. In such 
a case the residuals have a pattern on the sky indicative of 
the error made in the assumed principal-axis orientations. 

Given a reasonable starting point vo, the algorithm 
reliably detects non-axisymmetry in the "galactic plane." 
When the Sun-centre line does not lie too exactly in a princi- 
pal plane of the model, the shape of the recovered isodensity 
surface is very good away from the edge of the bounding box. 



When the Sun-centre line lies nearly in a principal plane, the 
detailed shape of the recovered isodensity surfaces is liable 
to be unrealistically angular, although their elongation is 
correctly recovered. 

The algorithm does not perform well when the data 
are very noisy. In general worthwhile results are obtained in 
regions for which the signal-to-noise is unity or greater. For- 
tunately, the quality of the results obtained at small radii, 
where the data may be expected to have high signal-to-noise, 
is not seriously undermined by the fact that every line of 
sight has to pass through low-density regions for which reli- 
able data will not be available. 

It will be interesting to see whether uncertainties due 
noise and dust absorption can be reduced to a level which 
allows this algorithm to be applied to the COBi?/DIRBE 
near-infrared data. 
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